Identification of potential modifier genes in Chinese patients with Wilson disease

Abstract The mutations in modifier genes may contribute to some inherited diseases including Wilson disease (WD). This study was designed to identify potential modifier genes that contribute to WD. A total of 10 WD patients with single or no heterozygous ATP7B mutations were recruited for whole-exome sequencing (WES). Five hundred and thirteen candidate genes, of which the genetic variants present in at least two patients, were identified. In order to clarify which proteins might be involved in copper transfer or metabolism processes, the isobaric tags for relative and absolute quantitation (iTRAQ) was performed to identify the differentially expressed proteins between normal and CuSO4-treated cell lines. Thirteen genes/proteins were identified by both WES and iTRAQ, indicating that disease-causing variants of these genes may actually contribute to the aberrant copper ion accumulation. Additionally, the c.86C > T (p.S29L) mutation in the SLC31A2 gene (coding CTR2) has a relative higher frequency in our cohort of WD patients (6/191) than reported (0.0024 in gnomAD database) in our healthy donors (0/109), and CTR2S29L leads to increased intracellular Cu concentration and Cu-induced apoptosis in cultured cell lines. In conclusion, the WES and iTRAQ approaches successfully identified several disease-causing variants in potential modifier genes that may be involved in the WD phenotype.


Introduction
Wilson disease (WD, OMIM entry: 277 900) is a rare autosomal recessive genetic disorder associated with various mutations in the ATP7B gene. It affects 1 in 30 000 individuals with a carrier frequency of 1 in every 90, the common symptoms including abdominal pain, jaundice, personality change, etc. The diagnostic algorithm is based on a diagnostic index ("Leipzig" score) proposed by an expert panel 1 . This score includes clinical findings (Kayser-Fleischer corneal ring), biochemical findings (low serum copper and ceruloplasmin concentrations, and increased urinary copper excretion), and molecular features (detection of biallelic ATP7B pathogenic variants). The mainstay therapy for WD is copper chelation therapy with penicillamine or trientine.
Previous research has shown the genotype of patients with WD varies between different populations. According to our previous research, in China, about 66.2% of patients with WD have compound heterozygous variants, and only 12.3% are homozygotes. 2 However, recent studies have indicated that approximately 12-20% of patients with WD have only one or no heterozygous mutation, whether from China or from Western countries. These patients also have some typical WD phenotypes such as increased urinary copper excretion and low serum ceruloplasmin concentration. 2,3 To date, molecular testing methods used to diagnose WD in the clinic have mainly relied on Sanger sequencing. 4,5 However, this approach is tedious and ineffective when sequencing large DNA fragments or searching for disease-causing variants in potential modifier genes. Whole-exome sequencing (WES) provides an efficient approach for identifying large fragment deletions and rare mutations 6,7 , and provides an accurate way to differentiate between systemic errors and actual biological variations. 8,9 The application of WES has already exhibited superiority to Sanger sequencing in the clinical genetic diagnosis of several kinds of inherited diseases. [10][11][12] A modifier gene is a gene that, when mutated, does not produce an effect on its own, but may produce or enhance clinical manifestations when coupled with another genetic mutation. 13 Several genes have been proposed as potential modifier genes for WD, including ATOX1 (antioxidant copper chaperone 1) 14 , COMMD1 (copper metabolism domain-containing 1) 15 , XIAP (X-linked inhibitor of apoptosis) 16 , and HFE (hemochromatosis gene). 17 However, none of these candidates has been confirmed by additional research.
Copper is an essential micronutrient required for oxygendependent enzymes. However, copper in excess is a noxious agent due to the oxidative capacity in biological systems. The opposing forces is balanced by chaperones and membrane transporters, which control copper uptake and moving copper across membranes. 18,19 The high-affinity Cu transporter, copper transport protein, CTR1, is responsible for the majority (∼70%) of Cu intake. 20 Upon entry, Cu may bind transiently to glutathione and be retrieved by Cu chaperones. Cytosolic Cu chaperone for SOD1 facilitates Cu incorporation, as well as the formation of a disulfide bond. 21 ATOX1 binds Cu(I) within the Cys-Gly-Gly-Cys metalbinding site, increasing the delivery of copper to the secretory pathway. 22 The P-type ATPase transporters, ATP7A and ATP7B, regulate cytoplasmic copper by pumping them out of cells or into the endomembrane system. In the trans-Golgi network (TGN), ATP7B transfers Cu to the ceruloplasmin, which is then secreted across the basolateral membrane. Cu elevation triggers trafficking of ATP7B from the TGN toward the apical membrane, thus facilitating export of excess Cu. 19 The CTR1 and CTR2 copper transporters are influx Cu(I) transporters in mammalian cells. CTR1 has a much higher affinity for copper than CTR2 and largely resides in the plasma membrane. CTR2 has a lower affinity for copper and is mostly associated with vesicular structures inside the cell. Only part of the CTR2 resides in the plasma membrane 23 , but it plays a role in stabilizing a truncated form of Ctr1 lacking the extracellular domain. Retention of this domain in mice or cells lacking Ctr2 enhances copper uptake. 24 Generally, CTR1 mobilizes Cu(I) from the vacuole when Cu(I) is scarce, but the function of CTR2 is far from fully understood. For the WD patients with typical clinical symptoms but only single or no heterozygous ATP7B mutation, it is plausible to hypothesize that the copper transporters (CTR1 and/or CTR2) might act as modifier genes to influence the WD phenotype.

Diagnostic criteria
A total of 201 patients with WD were recruited from the China Registry of Genetic/Metabolic Liver Diseases (CR-GMLD) Group. After laboratory and genetic tests, diagnosis of WD was established in accordance with the Leipzig score (≥4), and 10 patients with only one or no ATP7B mutation were selected for WES. An additional cohort of 191 patients with WD was recruited to determine the frequency of the CTR2 p.S29L mutation.

Sanger sequencing
Candidate pathogenic mutations in patients with WD and their immediate family members were identified by Sanger sequencing using an ABI 3730 analyser (Applied Biosystem). All of the 21 exons and at least 50 bp of intronic flanking sequences of the ATP7B gene were amplified, using primers reported previously. 2 Sites of variation were identified through comparing DNA sequences with the corresponding GenBank reference sequences (GRCh38, https://www.ncbi.nlm.nih.gov/genbank/) using Mutation Surveyor software.

DNA library preparation and sequencing
Peripheral blood (approximately 2 mL) was collected from the patients, and genomic DNA was extracted using the TIANamp Blood DNA Kit (Tiangen Biotech Co., China). DNA libraries were prepared according to Illumina's (USA) protocol. In brief, 3 μg of genomic DNA was fragmented by sonication (Covaris S2, USA). An "A" was then ligated to the 3 -end of the fragmented DNA. The sample was size-selected aiming to obtain segments of 150-∼250 bp for PCR amplification. Enriched libraries were then sequenced using an Illumina NextSeq 500 sequencer (Illumina) for paired-end reads of 150 bp.

Cell culture
HEK-293T, BEL-7402, and Huh7 cell lines were purchased from National Collection of Authenticated Cell Cultures (Shanghai, China), and incubated at 37°C with 5% CO 2 in a humidified atmosphere. The cell lines were maintained in 1640 (for BEL-7402) or Dulbecco's Modified Eagle Medium (for Human Embryonic Kidney Cells 293T and Huh7) medium (Life Technologies, USA) supplemented with 10% (v/v) fetal bovine serum (Life Technologies) and antibiotics (100 IU/mL penicillin and 100 μg/mL streptomycin).

Protein extraction and quantification
Huh7 cells (including ATP7B-knockdown and negative control cell lines) were incubated with 100 μM CuSO 4 for 0, 24, and 48 hrs. Proteins were extracted using lysis buffer (0.1 M Tris/HCl buffer containing 7 M urea, 2 M thiourea, 0.1% phenylmethylsulfonyl fluoride, and 65 mM dithiothreitol, pH 8.0), and sample concentration was determined using a Bradford protein assay kit (Thermo, USA); 50 μg of protein were used for each sample. The detailed protocol for isobaric tags for relative and absolute quantitation (iTRAQ) analysis has been published previously. 26

Liquid chromatography and tandem mass spectrometry analysis
Liquid chromatography and tandem mass spectrometry (LC-MS/MS) was performed by the National Center for Protein Sciences (Beijing). In brief, protein samples were precipitated with acetone-trichloroacetic acid and trypsin digested to generate proteolytic peptides. Proteolytic peptides were then labeled with iTRAQ reagents. The iTRAQ labels 113-118 were used to separately label the samples at various time points, and labeled samples were pooled before further analysis. The combined peptide mixtures were analysed using LC-MS/MS for both identification and quantification. The "abundances" indicated the quantification of the protein, and the fold change was calculated by comparing the abundances of different samples.

The wild-type and p.S29L (c.86C > T) mutant CTR2-expressed cell lines
Wild-type (WT) and mutant (Mut) human CTR2 cDNA sequences were chemically synthesized by Tsingke Biotechnology Co., Ltd (Beijing, China) and cloned into the pCDH-EF1-MCS-T2A-copGFP (pCDH) vector (MiaoLing Plasmid Platform, China); HEK-293T (National Collection of Authenticated Cell Cultures) lentivirus packaging cell lines were used to generate virus particles. Following co-transfection of pCDH, psPAX2, and pMD2.G plasmids (MiaoLing Plasmid Platform) in to 293T cells, replication-incompetent virions (CTR2 WT and CTR2 Mut ) were released into the culture medium supernatant, and the supernatant were collected after 48 and 72 hrs. The procedure was referred to the instructions provided by addgene (https://www.addgene.org/protocols/ lentivirus-production/). BEL-7402 and Huh7 cells (National Collection of Authenticated Cell Cultures) were co-incubated with CTR2 WT and CTR2 Mut lentivirus for 72 hrs. The copGFP + cells were selected using a flow cytometer, and were cultured as stable WT and Mut cell lines (BEL-7402 CTR2-WT/Mut and Huh7 CTR2-WT/Mut ).

Flow cytometry
Huh7 CTR2-Con/WT/Mut cells were seeded onto 6-well plates overnight. The cells were treated with 200 μM of CuSO 4 for 24 hrs, then 1mM D-PEN was added in the freshly changed culture medium for further incubation. After 24 hrs, the cells were collected and washed twice with PBS. The cells were stained with Annexin V-PE/7-AAD (PE Annexin V Apoptosis Detection Kit I, BD Biosciences, USA) on ice for 15 min, flow cytometric analysis was conducted, and apoptotic fractions were determined.

Quantification of copper ion concentration
Huh7 CTR2-Con/WT/Mut cells were seeded onto 100 mm cell culture dishes, and supplemented with 200 μM of CuSO 4 and/or 1mM D-PEN. After 24 hrs, the cells were harvested and washed twice by 2 mL PBS. The cell number was calculated by TC10 automated cell counter (Bio-Rad, USA). The quantification of copper ion concentration was measured using QuantiChrom TM Copper Assay Kit (BioAssay Systems, USA).

Statistical analysis
All in vitro experiments were repeated at least three times. Twotailed unpaired Student's t-tests and one-way ANOVA were used to evaluate differences between groups. GraphPad Prism 7 and MedCalc were used for statistical analyses. All quantitative data are presented as mean ± SD, and P < 0.05 was considered statistically significant.

Genetic and clinical profiles of the WD patients
Sanger sequencing was performed in our WD patient cohort recruited from the CR-GMLD Group. Single or no ATP7B mutations were detected in 10 of the patients. These patients were selected for further WES. The pathogenicity of the variants was analysed using prediction tools including Polyphen-2, MutationTaster, and SIFT (Table 1). Sanger sequencing identified seven different pathogenic variants in CDS regions, including the most common ATP7B variant p.R778L (3 in 10 patients), five other missense mutations (p.A874V, p.T935M, p.H1019P, p.V1106I, and p.N1270S), and one small In/Del mutation (p.F699-). In patient S3, no pathogenic variant was detected. SNPs were identified in patient S2 (p.S406A, p.V456L, p.K832R, and p.R952K), S3 (p.S406A, p.V456L, p.K832R, The clinical profiles of the 10 WD patients, including six males and four females from 10 unrelated families, are summarized in Table 2. In general, these patients had a high Leipzig score (≥4) and a low CP level (<0.2 g/L). Seven of 10 patients were hepatic phenotypes, including abnormal liver biochemistry and cirrhosis, and a hepatic and neurological mixed presentation was observed in the remaining three patients. In addition to the hepatic phenotypes, typical symptoms included tremors and involuntary movements.
The WES was performed by MyGenostics Inc. (Beijing, China). In addition of the pathogenic variants we detected by Sanger sequencing, in patient S5, we identified large fragment deletions (chr13:52 544 737-52 548 212 and chr13:52 549 161-52 585 503) in one of the ATP7B alleles. These mutations were identified by WES, but failed to be detected by Sanger sequencing (Table 1). Pedigree analysis revealed that the patient was a compound heterozygote.

Discovering potential modifier genes in nine WD patients
We analysed the WES data of the remaining nine patients with WD. Variants with a frequency <0.01 in 1000 Genome, ESP6500, ExAC_ALL, and ExAC_EAS databases were selected for further analysis. Multiple function prediction software, including Polyphen-2, MutationTaster, and SIFT, were used to predict the pathogenicity of nonsynonymous variants with uncertain clinical significance. Five hundred and thirteen candidate genes, of which the genetic variants present in at least two patients, were finally identified (Supplemental File S1). These candidate genes included 23 transcription factors and two splicing factors (Table 3), and 14 metal ion transport protein encoding genes [including SLC31A2 (encoding CTR2); Table 3]. These transporters might have synergistic or competitive effects on copper transport. However, no metal ion metabolism-related biological processes were significantly enriched (Supplementary Fig. S1). No pathogenicity variants were identified in those previously reported modifier genes, including ATOX1, COMMD1, XIAP, and HFE.
We next examined which proteins might be involved in copper transfer and metabolism processes, or affected by copper stimulation. The iTRAQ was performed to identify proteins differentially expressed between normal and Cu-stimulated cell lines. We generated the ATP7B-knockdown stable Huh7 cell lines by lentivirus.
The Huh7 ATP7B-NC and Huh7 ATP7B-KD cells were co-incubated with 200 μM of CuSO 4 for 0, 24, and 48 hrs, then the cells were harvested and protein samples were digested to generate proteolytic peptides. The proteolytic peptides were labeled with iTRAQ reagents. Protein expression patterns were quite similar in Huh7 ATP7B-NC and Huh7 ATP7B-KD cell lines ( Fig. 1A and B). Proteins with 1.2-fold increase or 0.83-fold decrease and P-value (<0.05) [27][28][29] were considered differentially expressed. Up-regulated (n = 363) and downregulated (n = 949) proteins were identified (Supplemental File S2). KEGG analysis revealed that the metabolism-related pathways were highly enriched by those down-regulated proteins; while the up-regulated proteins were mainly functioning in nervous system disease-related pathways, including Alzheimer disease, Huntington disease, and Parkinson disease-related pathways (Fig. 1C), which consistent with the typical clinical features of WD (liver abnormalities and neurological symptoms).
Thirteen genes/proteins were identified in both WES and iTRAQ data, including TTN, LAMA5, PRSS1, PLEKHA7, CDC27, DYNC2H1, FBLIM1, FDXR, MRTFB, NLRP13, PISD, NOP53, and PIH1D1 (Table 4). Only PLEKHA7 was reported to regulate the localization and function of ATP7A to promote copper extrusion in elevated copper through WW-PLEKHAs (PLEKHA5, PLEKHA6, PLEKHA7)-PDZD11 complexes 30,31 , no reports published to date indicate that the rest proteins are involved in copper ion intercellular transportation or metabolism processes. The molecular basis through which these proteins may be involved in WD requires further investigation. For some of the potential copper responsive proteins reported previously, the information from iTRAQ, WES, and literature search, were combined together in Table 5.

The CTR2 p.S29L variant promotes the Cu-induced apoptosis in BEL-7402 and Huh7 cells
CTR2 is an influx Cu transporter in mammalian cells. Because of the low affinity to copper ion, the function of CTR2 is still not fully understood. Our WES analysis identified a c.86C > T missense mutation in exon 3 of the SLC31A2 gene. An additional 191 patients with WD and 109 healthy donors were recruited for detecting the SLC31A2 c.86C > T (CTR2 p.S29L) variant using Sanger sequencing ( Fig. 2A). Our results showed a much higher SLC31A2 c.86C > T frequency in our WD cohort (6/191, 0.0314) than in East Asian (0.0024 in gnomAD database, P < 0.0001, test for one    Multiple function prediction software, including Polyphen-2, Mu-tationTaster, and SIFT, indicated that this single-base substitution mutation has strong pathogenicity ( Table 3).
We preliminarily assessed the pathogenicity of CTR2 p.S29L in cultured cell lines. Considering the low expression of endogenous CTR2 in BEL-7402 and Huh7 cells, we constructed BEL-7402 CTR2-WT/Mut and Huh7 CTR2-WT/Mut stable cell lines by directly overexpression of CTR2-WT/Mut genes using lentivirus (Fig. 2B). The cells were seeded onto 96-well plates and co-incubated with 200 μM of CuSO 4 for 24 hrs, then 1 mM D-PEN was added in the freshly changed culture medium for further incubation. CTR2 S29L caused a greater amount of Cu-induced apoptosis in both BEL-7402 CTR2-Mut and Huh7 CTR2-Mut cells after 24 hrs (Fig. 2C). The Cu-induced phenotype could be rescued by 1 mM of D-PEN at 48 and 72 hrs, in cells with negative control, WT, or mutant CTR2. The apoptosis of the Cu-stimulated cells was detected by flow cytometry. After co-incubation with 200 μM of CuSO 4 for 24 hrs, 19.935% of cells were apoptotic in Cu-stimulated Huh7 CTR2-Mut cells (15.643% in Huh7 CTR2-WT cells), and 1mM of D-PEN reduced the proportion to 12.730% (10.408% in Huh7 CTR2-WT cells) (Fig. 3). The increased Cu-induced apoptosis might be due to the enhanced uptake of the copper ion. We further assessed the intracellular Cu concentration in Huh7 cells, with CuSO 4 and/or D-PEN stimulation. CuSO 4 incubation caused an enhanced uptake of the copper ion in Huh7 cells, and D-PEN was able to decrease the intracellular copper concentration (Fig. 2D).
We further used siRNA pools to knock down ATP7B gene in Huh7 cells that express NC, WT, and mutant CTR2. The knockdown of ATP7B gene caused varying degrees of apoptosis in Cu-stimulated Huh7 cells, including the cells that with only endogenous CTR2. Although the cells with WT CTR2 also revealed enhanced apoptosis than negative controls, the cells with mutant CTR2 showed significant reduced survival rate compared with their WT counterparts at 24 hrs (Fig. 2E).

Discussion
Although with only single or no heterozygous mutations, some of the patients may still exhibit severe WD manifestations, they do not obviously differ from those with pathogenic homozygous or compound heterozygous mutations. In the current study, all of the 21 exons and at least 50 bp of intronic flanking sequences of the ATP7B gene were amplified by Sanger sequencing, and the multiplex ligation-dependent probe amplification test was performed to make sure that there was no large deletions in the promotor region of the ATP7B gene (data not shown). To our knowledge, no pathogenic mutations in deep intronic sequences that cause WD have been reported. Thus, the phenotypic variability and incomplete penetrance observed in WD suggests the involvement of modifier genes. To date, rare variants in several genes have been associated with WD, and a recent report using WES showed an increased and reduced risk of neurological presentation by variants in ESD and INO80 genes, respectively. 32 Furthermore, rare APOE and MBD6 variants are associated with lower risk of early onset WD. 32 Although the search for WD modifier genes is difficult, it is worth pursuing. The identification of these modifiers may provide a better understanding of the biological pathways that WD indirectly affects, and could also direct therapy possibilities for patients with WD. 33 Several strategies have been proved effective in identifying modifiers, including linkage analysis, association studies, biology-driven approaches, and systematic screens. 33 WES captures and sequences the entire exome; in the human genome, WES covers approximately 22 000 protein-coding genes. 34 This approach can also be used to capture low-frequency copy number variations or complicated In/Del mutations. 34 In this study, WES identified two large fragments deleted in the ATP7B gene in one of our patients, which was not detected by Sanger sequencing. However, a WES strategy might also be disappointing when searching for potential modifier genes, especially where limited samples are available. Here, we identified 513 potential modifiers with genetic variants in at least two patients, but no metal ion metabolismrelated biological processes were significantly enriched. The genes with potential mutations within the 3 or 5 UTRs (untranslated regions) which encode copper-binding proteins, may serve as modifier genes. The targeted next-generation sequencing enables rapid identification of genetic variations of a series of gene sets, both in exonic, intronic, and UTR regions. Thus, using the targeted next-generation sequencing for rare pathogenic mutations of copper metabolism-related genes will be helpful in searching the potential modifier genes of WD. Therefore, it may be wiser to focus on a more limited number of genes, the so-called candidate genes, on the basis of biological knowledge or molecular experiments.
To date, only a limited number of the variants identified in "modifier genes" have been confirmed by molecular experiments or other research. Using WES and iTRAQ approaches, we identified 13 candidate modifier genes, rare variants of which were supposed to be highly involved in the WD phenotype. However, the function annotations in the UniProt database (https://www.uniprot.org/) did not indicate a relationship between those candidates and WD phenotype. To our knowledge, so far, only PLEKHA7 was reported to regulate the localization and function of ATP7A to promote copper extrusion in elevated copper through WW-PLEKHAs (PLEKHA5, PLEKHA6, PLEKHA7)-PDZD11 complexes 30,31 ; no reports published to date indicate that the rest of the proteins are involved in copper ion intercellular transportation or metabolism processes. The expression levels of those proteins were changed under copper stimulation, but we are still not sure whether these candidates are eventually involved in the copper metabolism process, and if so, how they contributed to the WD phenotype was also largely unknown.
Although the differential expression of CTR2 was not detected in our iTRAQ analysis, the CTR2 p.S29L may also lead to the change of function of WT CTR2 even if no mutation-induced degradation happens. Although most of the CTR2 is localized in endosomes and lysosomes 23 , the biogenesis of the CTR1, the main Cu(I) influx transporter in mammalian cells, requires the existence of CTR2. 24 CTR2 might modulate Cu(I) uptake by controlling the expression or function of CTR1, and primarily functions as a regulator of influx, but has a limited effect on Cu efflux. 35,36 Here we reported a relative high patient frequency (6/191, 0.0314) of the CTR2 p.S29L (SLC31A2 c.86C > T) variant in our cohort of 191 WD patients, and even higher in atypical WD patients (2/9, 0.222), suggesting a specific effect of CTR2 p.S29L on atypical WD patients without the classical ATP7B biallelic mutation. To date, no reports have indicated a direct association between CTR2 and WD, but our study indicated that CTR2 p.S29L may be an important pathogenetic variant in Chinese WD patients, especially in those with atypical disease-causing genotype of ATP7B gene. However, because of the limited number of the patients carrying CTR2 p.S29L variant, we were not able to observe if they had any differences in clinical manifestations with other WD patients. In our study, when incubated with CuSO 4 solution, increased level of apoptosis was observed in CTR2 S29L -overexpressed cell lines. These results suggest that the CTR2 p.S29L variant might cause more copper accumulation than does the WT protein, and that CTR2 p.S29L may contribute to the clinical manifestations of WD. However, whether CTR1 is also involved in this process remains unclear, and more in-depth molecular research is urgently needed.